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X  ^  Abstract 


In  Part  II  of  this  paper »  the  approximation  theory  developed  in  Part  I 
“^til^is  used  to  determine  the  piecewise  polynomial  approximability  of 
^  solutions  of  elliptic  problems  on  polygonal  domains  in  Rt  and  polyhedra 
in  R'G.  From  these  estimates ,  convergence  orders  for  the  p-version  of  the 
finite  element  method  applied  to  such  problems  are  readily  obtained.  The 
critical  issue  is  the  approximation  of  the  singularities  which  occur  at  the 
ron-smooth  parts  of  tl^  domain  boundaries. |  It  is  seen  that  the  estimates  of 
''  tllj  involving  the  weighted  SuUulev  spacfes^  Zf  are  well-suited  for 

treating  such  singular  functions,  yielding  directly  the  optimal  approximation 
.degree,  up  to  an  arbitrarily  small  e. 

f'fcjmerical  results  for  two  problems  from  two-dimensional  linear  elasticity 
are  also  presented.  The  computations  show  that  the  predicted  order  of 
convergence  is  achieved  even  for  low  values  of  p.  Moreover,  in  contrast  to 
the  usual  h-version  of  the  finite  element  method,  the  point  at  which  the 
p-version  enters  the  asymptotic  range  does  not  depend  on  problem  parameters 
such  as  the  Poisson  ratio. 


(\ 

\ 

\ 


Accession 

ntis  gra&i 

DTIC  TAB 

Uri 'announced 


- - 

oy. - 

Distri 

Aval 

Dlst 

fi'i 

hilt  Ion/  _ i 

Lab 11 it 

Avail  * 
Spec! 

y  Codes _ 

lad/or 

.al 

_ 

1.  Introouction 

In  this  paper,  the  approximation  theory  developed  in  [11]  for  the 
p-version  of  the  finite  element  method  is  applied  to  elliptic  problems  in 
polygonal  and  polyhedral  domains.  Assuming  that  the  p-version  (see  e.g.  [6]) 
is  used  to  compute  an  approximate  piecewise  polynominal  solution  to  some 
elliptic  problem,  the  question  at  issue  is  the  rate  at  which  the  approximate 
solution  will  converge  to  the  true  solution  in  the  appropriate  energy  space  as 
the  polynomial  degree  p  is  increased.  It  is  well-known  that  if  a  coercivity 
or  "inf-sup"  condition  [2,10]  exists  for  the  given  elliptic  operator,  then 
this  question  reduces  to  a  purely  approx imation-theoretic  problem  of  the  kind 
treated  in  [11],  The  purpose  of  the  present  paper  is  to  show  how  the  results 
of  [11]  may  be  used  to  determine  the  approximability  of  solutions  of  some 
model  problems  in  the  usual  Sobolev  spaces  by  piecewise  polynomials  satisfying 
appropriate  boundary  and  conformality  conditions. 

As  with  any  approximation  result,  the  degree  of  approximation  for  the 
p-version  is  determined  by  the  regularity  of  the  function  being  approximated. 
The  regularity  of  solutions  of  elliptic  problems  on  smooth  domains  is 
classical  [1,19],  unfortunately  few  problems  of  practical  interest  fall  into 
this  category.  For  elliptic  problems  on  non-smooth  domains  such  as  polygons 
and  polyhedra,  one  must  turn  to  the  regularity  theory  developed  e.g.  in 
[12-18].  Essentially,  these  results  show  that  solutions  on  such  domains  may 
be  decomposed  into  the  sum  of  smooth  functions  and  functions  which  possess 
singular  derivatives  at  the  corners  (and  edges  in  R^)  of  the  domain.  It 
suffices  therefore  to  consider  the  approximability  of  these  singular 
functions,  whose  form  can  usually  be  given  very  explicitly.  It  has  been 
previously  noted  [6]  that,  unlike  the  h-version  estimates,  the  approximation 
results  for  such  singular  functions  as  obtained  from  estimates  involving  the 
usual  unweighted  Sobolev  spaces  are  not  optimal  for  the  p-version.  As 


will  be  seen  in  the  following  sections,  however,  the  singular  character  of 
these  functions  is  well-distinguished  by  the  weighted  Sobolev  spaces 
introduced  in  [11],  and  optimal  results  are  obtained  by  the  estimates  of 
[11].  This  is  clearly  due  to  the  close  relationship  between  the  weighted 
spaces  and  polynomial  approximation  indicated  in  [11]. 

In  Section  2  the  notation  and  main  approximation  results  of  [11]  are 
reviewed.  Section  3  addresses  the  approximation  of  the  solution  of  a  model 
elliptic  boundary-value  problem  in  a  polygon.  This  two-dimensional  result  has 
been  previously  obtained  in  [6]  by  a  different  method  not  involving  weighted 
spaces.  The  advantage  of  the  current  approach  lies  in  the  fact  that  more 
general  kinds  of  singularities  can  be  analyzed  using  the  weighted  spaces 
Z  than  with  the  techniques  of  [6].  In  particular,  one  may  allow  a 
continuum  of  singularities  along  part  of  the  bounoary  of  the  Domain.  Tnis 
situation  occurs  in  three-dimensional  problems  along  edges  of  polyheoral 
domains  as  described  in  Section  4,  where  some  new  approximation  results  are 
obtained.  In  Section  5,  some  numerical  results  are  presented  for  a  sample  of 
problems  in  two-dimensional  linear  elasticity.  These  were  computed  by  the 
research  program  COIET-X  [7]  developed  by  the  Center  for  Computational 
Mechanics  at  Washington  University  in  St.  Louis. 


2.  Review  of  the  Approximation  Results 

The  following  is  a  summary  of  the  notation  and  main  results  of  [11]. 

For  each  positive  integer  n,  let  I*^  =  {x  =  (x^,  ...,  x^^):  -1  <  <1, 

1  ^  i  _<  n  i.  For  each  non-negative  real  number  s  such  that  s  ^  +  an  integer , 


define 


Z®(l")  =  {u:  Hull  ^  „  <  •} 

Z®(l") 


where,  if  s  =  k  an  integer,  then 


and  if  s  =  k  +  B  with  k  an  integer  and  0  <  B  <  1  (B  ^  1/2),  then 


I 

i=l  ^n-l 


,,  ^.2x5/2  3  u  . 

(1  —  t  )  ^  •••>  ^i-1  ’ 

^i 


-  (1  -  (X  X  T  X  X  ) 


•  It  -  dt  dxdx,  ...dx.  ,  dx.  .  ...dx 

1  1-1  1+1  n  I 

Let  a  be  a  Domain  in  such  that  there  exists  a  triangulation  A  of 
ft  into  open  n-sitnplices  ft^^,  i=l,  ...,  M.  Let  o\  v  =1,  ...,  N, 
oenote  the  vertices  of  A.  Since  it  is  convenient  to  be  able  to  refer  to  the 
vertices  of  a  particular  simplex  ft.,  for  i=l,  m  let  a.  1  £  j 
<,n  +  1,  denote  the  vertices  of  ft^.  Hence,  if  a  vertex  of  A  is 
also  a  vertex  of  ft..,  then  o'*  =  o.  .  for  some  j.  Let  S.  = 

i  1,J  D 

1  <  i  <M,  1  1  j  in  +  1). 

Consider  a  simplex  ft.  e  A  and  one  of  its  vertices  o.  . . 

^  ^  >  J 

Letting  e^,  ...»  ©p  be  the  vectors  based  at  j  and  terminating  at 
the  other  n  vertices  of  ft^,  define  in  the  parallelepiped 


n 

ft.  .  =  {  X  e  FT:  X  =  o.  .  +  Z  t  e„,  0  <  t  < 

mm  m 


1}  . 


Clearly,  ft.  c  ft  for  all  j=l,  ...,  n+1.  For  each  (i,j)  c  S,  choose  an 

1  ifj  a 

affine  mapping  T.  .  in  R^  which  maps  ft..  .  onto  and  o.  .  onto  the  point 
^fj  J  ^  fJ 


N  • 

Let  be  a  smooth  partition  of  unity  on  Si  such  that  for  each  \tel, 

...»  N,  supp  contains  the  vertex  and  supp  ri^  intersects  only  those 

closed  simplices  Si^  which  have  o^  as  a  vertex.  For  any  u  e  L2(Si)  and  (i,j) 


e  S^,  define  in  Si^  ^ 


i»J 


un  in  Q.  where  v  is  such  that  o.  .  =  o' 

V  1.  1 » J 

0  in  Si.  ./Si.  . 

*»J  i 


By  the  assumptions  on  rv^,  one  observes  that  u.  .  =  u  in  a  neighborhood 

^  r  I J 

of  o.  ,  and  u,  ^  =  0  outside  of  a  neighborhood  G  of  a.  .  such  that 

•••ij  •*-ij  iij 

G  ih  Si.  ./  s<.  = 

For  each  real  s  ^0  such  that  s  /  ^  +  an  integer,  let 

Z®(Si;A)  =  {u  :  (lull  <  ««>} 

Z®(Si;A) 

where 


lUH 


Z®(Si;A) 


•( 


(i,j)6S^ 


10.  .0 


rT^  .11^ 

i,J  zS(in) 


y/2 


For  each  non-negative  integer  p,  let  1  1  i  1  denote  the 

space  of  all  polynomials  on  Si^  of  degree  at  most  p,  and  let 
Pp(4i;A)  =  iu  :  ulj^  ePp(Si^),  1  <_i  <M} 

For  all  non-negative  integers  1  and  p,  set 
Z®  =  Z®(Si;A)  n  c^ri)  , 

Pp  (Si;A)  =  Pp(Si;A)  n  c^H)  , 

Jt  — 

where  C  (G)  denotes  the  set  of  all  functions  which  along  with  their 
first  H  derivatives  are  continuous  on  Si.  Pp(Si;  A)  is  then  the 
set  of  all  functions  in  Pp(Si;A)  which  along  with  their  first  ^ 
derivatives  are  continuous  across  the  common  boundaries  of  adjacent  simplices 


’SiSi 


of  A.  It  is  shown  in  [11]  that  if  0  <_  ii.  <  then  the  analogous 
statement  holds  for  2^(il;A)  and  Z^(tt;A). 

Recalling  that  il  is  the  union  of  n-simplices  fl. ,  let  r  ,  k=l, 

K,  oenote  the  faces  of  fl,  i.e.,  the  (n-l)-simplices  whose  union  is  the 
boundary  of  ii.  Fix  s  ^0  and  let  £?»  be  the  largest  integer  strictly 
less  than  .  For  each  k=l,  K,  let  be  a  subset  of  {0,  1, 

10^},  Denoting  by  D  the  family  {Dj^;  1  <k  _<  K),  define 

Z®(D;A,D)  =  {u  e  Z^(a;A)  :  for  k=l,  K, 


— -  =  0  on  E  for  each  m  e  D,  } 
3n^  k  k 


Pp(fl;A,D)  =  {u  ePp(D;A)  :  for  k=l,  K, 


=  0  on  r.  for  each  m  e  D.  } 
^m  k  k 


where  — pr  =  u  and,  for  m  >  1,  — ^  is  the  order  derivative  of  u  with 

'<rP  “  an*" 

respect  to  n,  the  outward  pointing  unit  vector  normal  to  the  boundary  of  D. 

Letting  l-l®(fl)  be  the  usual  Sobolev  space  of  order  s  on  P  ana 
setting  Z®^(a;A,D)  =  Z®(fl;A),  the  main  result  of  [11]  may  now 
be  stated: 


Theorem  2.1  Let  s  and  s'  be  such  that  s  >  2s'  ^0  and  s,  2s'  ^  +  an 

integer.  Let  denote  the  largest  integer  strictly  less  than  If  u 

e  Z^(ft;A,D)  for  some  integer  I  satisfying  s'  -  |-  <  £.  <  £*  and  for  some  D,  then  for 

£ 

each  non-negative  integer  p  there  exists  cpp  e  Pp  (fl;A,D)  such  that  for  arbitrarily 
small  e  >  0, 

(2.1)  lU  -  (p  II  <  Cp’  ® 


k*3 H  yim  I.IHI.H l«  W«l''  I* «IV’ 


where  C  =  C(s,  s',  e)  is  independent  of  u  and  p.  Moreover,  if  u 
e  2^*(ft;A,D)  then  for  any  integer  t  >  s'  -  3/2  and  for  any 
non-negative  integer  p  there  exists  q)p  e  Pp(ft;Zi,D)  such  that 
(2.1)  holds. 


For  any  real  s  ^0,  define  D  =  {Dj^}  as  above  and  let 
H®(fi;D)  =  (u  e  :  for  k=l,  ...,  K, 

3f'\j 

=  0  on  for  each  m  c  ) 

3n 

The  following  [11, Theorem  3. A]  is  a  consequence  of  Theorem  2.1; 


Theorem  2.2  Let  s  and  s'  be  real  numbers  such  that  s  >  s'  ^0.  If 
u  e  H^(J2;D),  then  for  each  non-negative  integer  p  and  each  integer  £  > 
s'  -  3/2  there  exists  <Pp  e  Pp(ft;A,0)  such  that  for 
arbitrarily  small  e  >  0, 


-  9n»  s',  . 

P  («) 


-  s  +  s'  +  e 


Hull 


where  C  =  C(s,s',e)  is  independent  of  u  and  p. 


Remarks ; 

1.  Although  boundary  conditions  were  not  explicitly  included  in  [11,  Theorems 
3.3  and  3. A],  their  proofs  are  easily  modified  as  inoicated  in  [11]  to  obtain 
Theorems  2.1  and  2.2. 

2.  Theorem  2.2  was  previously  obtained  in  [6]  for  the  case  s'=l  and 

ul  /  rt\  ....  1 


iS 


iroximation  of  Solutions  of  Elliptic  Boundary-Value  Problems  in 


Polygonal  Domains 


Let  u  denote  a  polygonal  domain  in  R  with  bounaary  r.  Then  r  is 
the  union  of  a  finite  number  of  line  segments  1  £  <  L,  which  one 

may  assume  to  be  labelled  in  consecutive  order  following  the  positive 
orientation  of  F.  Assume  also  that  the  set  {1,  2,  ...,L}is 
partitioned  into  two  subsets  D  and  N,  and  let  F^  =  Fj^  and  Fj^  =  F^^. 

For  each  £,  1  £  fi-  £L,  let  n^  denote  the  outward  pointing  unit  normal 
to  Fjj^,  let  ^  if  i£D,  and  let  =  0  if  ieN. 

Consider  the  following  model  problem.  Given  k  a  non-negative  integer,  f 
e  hl'cn),  and  g^  e  ^*1/2  (r 

if  ieN,  find  u  such  that 


(3.1) 


-  A  U  +  U  =  f 


u  =  0 


an«  “ 


on  Fj,  ilcM 


Letting  H^(Ji;D)  =  {u  e  u  =  0  on  F^^},  (3.1)  is 

equivalent  (see  e.g.  [13])  to  the  following  variational  problem:  Find  u  e 

H^(J1;D)  such  that 

(3.2)  J  (^  •  W  +  uv)dx  -  Z  J  gj^  V  dF  =  J*  f  v  dx 


for  all  V  e  H‘^(r<;D) 


It  is  known  [13]  that  (3.2)  admits  a  unique  solution  u. 

The  present  objective  is  to  determine  the  approximability  of  the  solution 
u  in  H^(ft)  by  elements  of  Pp(f<;A,D)  where  A  is  some  triangulation  of  C.  In 


order  to  use  the  results  of  the  previous  section,  one  must  know  the  regularity 
of  u.  Since  f  e  it  might  be  expected  that  u  e 

However,  it  is  well-established  that  this  is  not  true,  in  general,  oue  to 
possible  singularities  at  the  corners  of  n  and  at  points  of  F  where  a 
change  of  boundary  condition  occurs.  Nevertheless,  under  certain  conditions, 
the  solution  u  may  be  decomposed  into  the  sum  of  a  function  in  and 

functions  which  possess  singular  derivatives  at  the  corners  and  points  of 
changing  boundary  conditions.  In  order  to  give  a  precise  statement  of  this 
result,  for  each  fi.,  1  £  £L,  let  o.  denote  the  common  endpoint  of  F.  and 

Fjul  ^1  “£>  ^  06'^ote  the  measure  of  the 

interior  angle  of  ii  at  a..  It  will  be  assumed  that  if  u).  =  n  then 
either  UdD,  Ju-l  eN  or  else  jleN,  il+l  eD  (with  the  obvious  modification  for 
i=L);  i.e.,  if  tt  then  the  boundary  conditions  in  (3.1)  change  from 
Dirichlet  to  Neumann  or  vice-versa  at  For  each  i,  1  1  t  ft,  let  (r^^, 
denote  polar  coordinates  based  at  such  that  6^=0  contains  the 
segment  Finally,  for  each  £,,  1  £  t  £L,  and  each  integer  m,  let 

6«*1  -  ’ 


The  following  is  a  consequence  of  [13]: 


k  k+1/2 

Lemma  3.1  Let  k  be  a  non-negative  integer,  let  f  e  H  (Ti),  and  let  o^e  h 

-  Bo  -  (k+l)u)„ 

(Fj^)  if  £.€N.  If - is  not  an  integer  for  any  Jt,  then  there 

k+2 

exists  w  e  H  (W)  and  numbers  c„  such  that 

x.,m 


£,  m 


^£,  m  ^  integer 


»i 

V, 


V 

■5 

■ » 


(3.3) 

^  0<X  ^  <k+l 
Jl,in 

^i,  m  “  integer 

is  the  unique  solution  of  (3.1),  where  the  are  smooth  functions  equal 
to  unity  near  and  zero  away  from  o^.  Moreover, 


(3.4)  iiwli 


‘  k+2  ^  ^  “^jL  m'  k  *  ^  k+1/9 

Hk+2(^)  x,,ni  £^k+l/2^P 


)  . 


) 


The  importance  of  Lemma  3.1  lies  in  the  fact  that  the  piecewise 
polynomial  approximability  of  u  may  now  be  determined  by  investigating  the 
approximability  of  each  term  in  the  right-hand  sioe  of  (3.3)  separately  (with 
a  slight  modification)  using  the  results  of  Section  2.  The  following  lemma  is 
required: 


Lemma  3.2  Let  J!.  be  such  that  1  £  <,  L,  let  e  c“(fO  be  such 

that  supp  r\^  is  contained  in  a  neighborhood  of  and  let 
00 

^  e  C  ([0,  w^]).  If  p  is  any  non-negative  real  number  and  if 
then  vTijj^  e  ZQ^^“^(ft;A)  for  arbitrarily  small  e  >  0. 

Pf :  Let  (r,6)  denote  a  polar  coordinate  system  centered  at  the  point  (1,1); 

2  2  1/2 

in  particular,  let  r  =  [(x^^-l)  +(x2-l)  ]  and  let  6  be  the  angle 
measured  from  the  semi-axis  {x:  x^  <  1,  X2  =  1).  It  suffices  to 
Show  that  if  ^  e  c“’([0,  ^]),  then  the  function 

v(r,e)  =  r^(Sn  r  ^^(0)  +  ^{Q))  ,  (r,  6)  e  I^  , 


N 
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Delongs  to  for  arbitrarily  small  e  >  0. 

To  this  end,  let  x  c  C*([0,“))  be  such  that 
0  for  0  <  t  <  1/2 

x(t)  = 

1  for  1  £  t  <  ® 

ana  for  each  6  >  0  let  X^(t)  =  Setting  =  vx^(r),  it  follows  by 

Leibniz*  rule  that  for  each  integer  k  >  2p  +  2  and  for  i  =  1,2, 


/  I  — I  U  -  x?)^  dx  <  c  if  I  I  I — i7-^  I  (1  -  X,)“  Ox 

■(2  ax^  ^  "  m=0-{2  ax”  " 

4-2k,2m  ^2u-2«^M-c 

V  m=0  ■'S/2  ■'S/2 


<C  . 


Furthermore,  since 


livJi  o  <  C  iivil  ^  , 

one  obtains  that  for  any  integer  k  >  2m  +  2, 

is  1  £ 

(3.6)  llvjl  .  ,  <C6^"2'^-^’2 

6  2k(j2j 


r  iv  -  v,i^  dx  =  r  iv(i  -  Xx) ^  dx 


.^1  w;  i 


For  each  t  >  0,  let 

K(v,t)  =  inf  (llv,  II  o  +  t  9  ) 

v=Vi+V2  L2(r)  ^  r(i  ) 

9/k 

If  0  <  t  <1,  then  by  choosing  6  =  t  it  follows  from  (5.6)  and 

(3.7)  that  for  k  >  2y  +2, 

K(v,t)  <  llv  -  V  .11  o  *  t  llv  Ji  .  „ 

‘  LjCI^)  ‘  Z^l^) 

2pf2-£ 

<  C  t  k 


Furthermore ,  for  t  1 , 


K(v,t)  <  iivli  „ 
U(I^) 


Hence,  for  arbitrarily  small  e  >  0, 


1  ¥-1 


St))^  ^  t  dt 


2  r 

Lo(I^)  *^1 


_  2vn-2-e 


dt  <  * 


This  implies  (see  e.g.  [9])  that  v  e  (L2(I^),  Z*^(I^))  2y+2-e 

k  ’  ^ 

and  the  result  follows  from  [11,  Theorem  2.1]. 

k  k+1/2 

Theorem  3.1  Let  k  be  a  non-negative  integer,  let  f  e  H  (ii),  let  Qn  ^  ^ 

ej^l-B,-{k+l)u), 

(Tjj^)  if  £eN,  and  assume  that - - - - is  not  an  integer  for  any  £. 

If  u  is  the  solution  of  (3.1)  then  for  each  non-negative  integer  p  there  exists 
e  Pp(I<;A,D)  such  that 

-  min{2X,  ,,  ...,  2X  ,,k+l}  +  e 

(3.8)  lU  -  \t<„ll  1  <Cp  (iifii  . 

P  H'^(fl)  H'^(fi) 


*  A. 


-  *  '  "".y. .  r*.  •  ^  - 


for  arbitrarily  small  c  >0,  where  C  =  C(e,k)  is  independent  of  u  and  p. 


V 


K'i 

w 

I 


l6 


Pf;  Let  w,  and  X£  oe  as  in  Lemma  3.1.  For  each  I 

and  m  such  that  0  <  <  k+1  and  ^  an  integer,  let 

^l,m  V  *  ■  ®Afl^  ^ 

For  each  £  and  m  such  that  0  <  X,  ^  <  k+1  and  X.  „  =  an  integer, 

£,m  £,m  ®  ’ 

define 


“‘tXTunrl 


(-1)"**^  u  (sin  6j^)  if  =  n,  £fl  e  N,  and  £eD 


(_1)"**^  «  (-cos  6jj^)  if  =  n,  id),  and  £+1  eN 


ano  set 


2.  Jr„e,)  «  r  »  {In  r.  cos(X,  „  e,  -  g,  ,)  -  e„  sin(X„  6  ,  -  J 


£,m'*  £’’'£'  *£ 


l,m  ''£  ^£+1'  m''  £  *"£+!' 

*  ^£,m^®£^^  ^£^^£’®£^  • 


Finally,  let 


(3.9) 


W  =  w-  £  ^  a  t  '  ^PtnXo 

^£  m  ■  integer 


0<X,  4<+l 

iuy  fT) 


2X,  +2-e 

£,tn 


Moreover  w  eH  (R)  o  h  (ft;D)  and  by  Lemma  3.2,  ^  ^  ^  ^  ’  (ti;A,D)  for 

arbitrarily  small  e  >  0.  Applying  Theorem  2.2  to  w  ana  Tneorem  2.1  to  tne 
z#  «.»  it  follows  that  for  each  non-negative  integer  p  tnere  exists  iL 

X|m  p 

cP°(m;o,D)  and  c  P?(f*;A,D),  0  <  X.  <  k  +  1,  such  that 

P  P  P  I»i  m 


ihV iiJ^Wid 


(3.10)  llw  -  “5, II  1  <  ilwii  .  ,  ,  and 

P  H^(a)  “  KT'^^Cfl) 


i(  iiiti)  I 


“2X«  +e 
£,111 


(3.11)  llz,  „  -  ilA‘'’'"^ll  <  Cp  *’"■  llz,  II  ,,  ,  ,  0  <  X.  <  k+1, 

where  C  is  independent  of  p.  Setting 


(3.12)  4.  .  ♦  .  £  ♦ 

P  P  0<X„  <k.l  P 

x,m 


(i,m) 


(3.8)  follows  from  (3.9)  through  (3.12)  together  with  (3.4)  and  the  triangle 
inequality . 


Remarks  1.  Theorem  3.1  for  the  case  Fq  =  r,  Fj^  =  0  is  proven 
in  [6,  Theorem  4.3].  The  techniques  employed  there,  however,  do  not  use 
weighted  spaces  to  obtain  the  result;  instead,  the  singular  functions  are 
approximated  directly. 

2.  Theorem  3.1  shows  the  doubled  order  of  convergence  obtained  by  using  the 
p-version  of  the  finite  element  method  instead  of  the  h-versioh  with 
quasi-uniform  refinement.  Letting  N  denote  the  dimension  of  Pp(fl;A,D) 
and  letting  h  denote  the  maximum  diameter  of  any  triangle  in  A,  it  holds  that 


(3.13) 


for  the  h-version  (quasi-uhiform  refihement), 
for  the  p-version. 


It  can  be  shown  (see  e.g.  [4])  that  functions  of  the  form  (3.5)  belong  to  the 
space  H®(ii()  iff  1  ^  s  <  X  +  1.  Letting 

M  =  min  1^^^ 

0<Xe  k+1 
£,m 


and  assuming  that  k  is  large  (which  it  usually  is  in  practice),  it  follows 
from  the  standard  h-version  approximation  theory  [2,10]  that  the  approximation 
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error  in  for  the  h-version  is  of  order  'w  n"  As  for 

the  p-version,  Theorem  3.1  together  with  (3.13)  yield  an  approximation  error 

of  oroer  p  for  arbitrarily  small  e  >  0. 


Results  such  as  Theorem  3.1  may  similarly  be  obtained  for  elliptic 
operators  other  than  that  of  (3.1)  provided  that  results  such  as  Lemma  3.1  are 
available  for  such  operators.  Consider  for  example  the  biharmonic  problem; 
given  f  e  H^(Q),  k  ^0,  find  u  such  that 


(3.14) 


^u  =  f  in  n 

u  =  3u/3n  =  0  on  r. 


Letting  H^(fl;D)  =  {  u  e  H^(ft)  :  u  =  au/an  =  0  on  r},  one 

2  1 

seeks  to  approximate  u  in  H  (ft)  by  elements  of  Pp(ft;A,D)  = 
F^(ft;D)  n  Pp(ft;A).  In  [183  it  is  shown  that  the  solution  u  of 
(3.14)  may  be  written  as 


u  =  w  + 


0  <  Re  <  a(k) 


Qi,m(rjl^X£(ri,ei) 


where 


(i)  w  e  H‘^'"^(ft;D), 


(ii)  the  X^  are  the  solutions  of  the  nonlinear  equation 


(iii)  the  functions  are  of  the  form 

=  £  Cmj(6i)  rjt  Sn 

for  some  smooth  functions  c^^^  j,  non-negative  integers  j,  and  real 
numbers  j  ^0,  and 

(iv)  is  a  smooth  cutoff  function  equal  to  unity  near  and 
vanishing  away  from  a^. 

Then  the  following  may  be  obtained  in  a  manner  similar  to  Theorem  3.1: 

Theorem  3.2  Let  f  e  H^(a)  for  some  non-negative  integer  k  and  let 
u  be  the  solution  of  (3.14).  Then  for  each  non-negative  integer  p  there 
exists  4ip  ePp(a;A,D)  such  that 

II  u  -  !|»p  II  H2(ft)  <  Cp“  2  min  ReXjj.^m  -  2,  k  +  2}  +  e  nfn 

SLf  m 

for  arbitrarily  small  e  >  0,  where  C  =  C(e,k)  is  independent  of  u 
and  p. 


4.  The  Approximation  of  Solutions  of  Elliptic  Boundary-Value  Problems  in 
Polyhedral  Domains 

'  3 

Let  id  denote  a  polyhedral  domain  in  R  with  boundary  P  and 


consider  the  following  model  problem;  Given  f  in  C  (fl),  find  u  such  that 


(4.1) 


■4  U  +  U  =  f 
u  =  0 


in  id, 
on  r  . 


Letting 


H^(C;0)  =  {u  c  u  =  0  on  r}  , 


fe 


w 

ck'M- 


it  is  well-known  that  (4.1)  admits  a  unique  variational  solution  u  e 

H^(Q;0).  In  order  to  determine  the  approximaOility  of  u  in  by 

elements  of  Pp(i2;A,D)  for  some  triangulation  4  of  ft,  the 

regularity  of  u  must  again  be  established.  By  classical  elliptic  regularity 

theory  [1,19],  the  solution  u  of  (4.1)  is  infinitely  differentiable  in  the 

interior  of  ft  and  up  to  the  smooth  faces  of  T.  However,  in  analogy  with 

the  previous  two-dimensional  case,  u  will,  in  general,  possess  singular 

derivatives  at  the  corners  and  edges  of  ft. 

Consider  first  the  behavior  of  u  near  an  edge  of  ft.  Let  y  denote  an 

open  segment  such  that  y  is  contained  in  the  interior  of  some  edge  of  ft. 

Letting  x  =  (Xj^,  X2,  x^)  denote  Cartesian  coordinates  in  R^,  assume 

without  loss  of  generality  that  y  is  contained  in  the  positive  x^-axis. 

Let  0)  denote  the  measure  of  the  interior  dihedral  angle  formed  by  the  two 

faces  and  r2  of  ft  which  intersect  along  y.  Let  (r,  <p,  Xj) 

denote  cylindrical  coordinates  about  the  x^-axis  such  that  r  =  (x^  + 

2  1/2 

X2)  and  such  that  <p  =  0  (resp.  =  u)  contains  F^^  (resp. 

2  — 

F2).  Denote  by  G  an  open  subset  of  R  such  that  G  x  y  c  ft  and  G  x 
y  neither  contains  any  corners  of  ft  nor  meets  any  eoges  of  ft  except  for 
the  edge  containing  y.  For  each  integer  m  >  0  and  each  integer  n  >  0, 
let 


mfl 

“  *  Sine  .) 


if  is  not  an  integer,  and  let 


mil  - 
+  2n 


m,n 


[in  r  sin(—  (p)  +  cos^—  (p)j 

U)  U) 


if  ^  is  an  integer.  The  following  result  is  contained  in  [16]: 
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Lewna  4.1  Let  f  belong  to  C*(n).  Then  for  any  positive  integer  k 

there  exists^  h  e  C*(Yi  and  e  C*(y)  such 

that 


14.4) 


u  = 


0  <  —  +  2n  <  k+1 

(I) 


Cm  „  2„  *  h 

in,n  ni,n 


in  G  X  Y* 


Remark  If  one  were  only  to  assume  that  f  e  H^(G),  Lemma  4.1  still  holds 

2  k+2 

except  that  in  this  case,  one  may  only  conclude  that  h  e  L  (y,  H  (G))  and 
k+1  -  M  -  2n  ^ 


c„  „  e  H  (y)  (see  [16]  and  the  references  contained  therein). 

m,n 

Od  V 

The  assumption  that  f  e  C  (G),  which  usually  holds  in  practice,  has 
been  made  in  order  to  help  simplify  the  presentation. 

Together  with  Lemma  4.1,  the  following  enables  one  to  apply  the 
approximation  results  of  Section  2  to  the  solution  u  of  (4.1)  in  a 
neighborhood  of  the  eoge  y. 


Lemma  4.2  Let  ^  e  C*([0,  w]),  let  n  e  C^^tG)  be 

such  that  supp  ti  c  G  x  y>  and  let  c  e  C*(y).  Let  v  be  a  function 

on  G  such  that,  in  G  x  y» 


C*(X;Y)  is  the  space  of  all  infinitely  differentiable  functions  on  X 
whose  values  lie  in  Y. 
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L2(X;Y)  is  the  space  of  all  functions  h  on  X  with  values  in  Y  such  that 
iih(x)iij  dx  <  •. 


v(r,q),Xj)  =  c(Xj)r^  (in  r  +  ^(9)) 

for  some  non-negative  real  number  y.  Then  vn  e  A,D)  for 

arbitrarily  small  e  >  0. 

Pf:  Since  v  is  smooth  in  the  variable  x^,  the  result  follows  from  Lemma  3.2. 


Consider  next  the  behavior  of  the  solution  u  of  (4.1)  in  a  neighborhood  of  a 
corner  of  U,  which  may  be  taken  to  be  the  origin  without  loss  of  generality.  Let 
Y  denote  an  edge  of  ft  which  has  the  origin  as  an  endpoint;  again  assume  that  y 
is  contained  in  the  positive  x^-axis  and  that  to  is  the  measure  of  the  interior 
dihedral  angle  of  ft  at  y.  Let  K  denote  the  intersection  of  ft  and  a  solid 
circular  cone  whose  vertex  is  the  origin,  whose  axis  contains  y,  and  whose 
half-angle  is  sufficiently  small  that  K  does  not  meet  any  of  the  edges  of  ft  which 
have  the  origin  as  an  endpoint  except  for  the  edge  containing  y.  Finally,  let 
denote  the  intersection  of  K  and  a  small  sphere  of  radius  about  the  origin. 

The  present  goal  is  to  understand  the  behavior  of  u  in  K^. 

To  this  end,  let  G  denote  the  section  of  the  unit  sphere  intersected  by  the 
infinitely  extended  cone  coinciding  with  ft  in  a  neighborhooo  of  the  origin,  and  let 
denote  the  intersection  of  K  and  the  unit  sphere.  Denote  by  0  <  ^ 

^2  i  • • •  the  eigenvalues  of  the  Laplace-Beltrami  operator  on  G  with  homogeneous 
Dirichlet  boundary  data  on  3G.  Let  (0,6,9)  be  the  usual  spherical  coordinate 
system  in  R^.  Applying  the  previous  two-dimensional  results  to  the 
Laplace-Beltrami  operator  on  G,  let  (6,9)  denote  the  singular  functions 
analogous  to  (4.2),  (4.3)  where  the  singularity  occurs  at  the  corner  6=0  of  G.  The 
following  is  proved  in  [16]: 


Lemma  4.3  Let  f  Delong  to  C  (U)  and  let  u  oe  the  solution  of  (4.1). 
Tnen  for  any  integers  k  and  k'  satisfying  0  <  k'  £  k,  there  exist  numbers 
Cm^n,i»  '^i  functions  g^,  h  such  that,  in  K^, 


u(p,  6,(|))  =  Z 


r  2 
L  1 


(4.5) 


c  .  p 
m,n,i  ^ 


- 1  I 


1  R  1 

-  T7  +-\,  A.+  -r 


^m,n^P^m,n^®’'P^  ^  ^  ^  "  %.(e,cp)  .  h(6,cp) 


■  2  \V  4 


where  g^  e  (G^)  and  for  arbitrarily  small  e  >  0, 

-k  +  j-  j  +  egj 

^  '™.n  ^ 


,  0  <  j  <  k'  , 


-  k  +  j  -  2  +  e  //r,  V  -  k'  +  2  _ 

p  h  e  L2  ((0,“);  H  (G^),  0  <  J  <k'  . 

The  following  lemma,  together  with  Lemma  4.3,  will  allow  the  application 
of  the  approximation  results  of  Section  2  to  the  solution  u  of  (4.1)  in  K^. 


Lemma  4.4  Let  rij^  e  C*(G)  be  such  that  supp  c  G^  and  let 
€  C*([0,<»))  be  such  that  supp  c  {p:  0  £  p  _<  p^}. 

(i)  If  V  is  a  function  on  S2  such  that,  in  K^, 

+  -Vu+So”^ 

(4.6)  p  ^  Sin  ^  £I-2((0,»);  L2(G))  , 

o<i6i  =  e^  +  ^  +  e^<v  , 

for  some  positive  real  numbers  y,,  yu  and  for  some  positive  integer  v,  then 
2  min(iJ,,W2)  +  2-c 

'^'^1’^  €  Zg  (ft;A,D)  for  arbitrarily  small  e  >  0. 

(ii)  If  V  is  a  function  on  ft  such  that,  in  Kq, 


f 


p 
I 

^  r 


If* 

|r 

m 


* 


■  :‘‘K 

-ifi.'.j 

mj 


n 

• '  '^i 


i  ^.-;  j,. 


4^ 

V/’fe 


M 


'f\, 

■■Xi 


^7^ 

iS'n 


’!»« 


-  ji  +  Bl  glii^ 


ap^aracp-^ 


€  L2((0,“);  L2(G)) 


01lll  =  6^  +  ^  +  &5lv 


2lh*2—£ 

for  some  positive  real  y  and  some  positive  integer  v,  then  vr^r^  ^  ^ 
(^A,D)  for  arbitrarily  small  c  >  0. 

Pf:  Assume  that  the  Cartesian  coordinate  system  (x^,  X2,  x^)  has  been 
transformed  by  an  affine  mapping  in  such  a  way  that  p  now  represents  the 
distance  to  the  point  (1,1,1)  and  that  6  is  the  angle  measured  from  the 
semi-axis  {x:  x^^  =  X2  =  1,  <  1).  To  prove  (i),  it  suffices  to 

show  that  a  function  v  in  satisfying  (4.6)  (with  the  new  definition  of 
(p,  8,9)  and  corresponding  redefinition  of  belongs  to 
2  min  (y,  ,1X5)  +  2-e  , 

Z  ^  d'^)  for  arbitrarily  small  e  >  0. 

To  this  end,  let  x  e  C*([0,*))  be  such  that 


X(t)  = 


0  for  0  <  t  <  1/2 
1  for  1  <  t  <  ~ 


and  for  each  6  >  0  let  X^Ct)  =  X^'^)*  Set  v^  = 
V  Xj(p)  Xj(P  sin  8).  Then,  for  i  =  1,  2,  5, 


,aL 

1  ‘•‘I  *i 


<  c  Z  6°^“' 


a’^xjip)  9^Xj(p  sip 


af^v,  2 

I  o  I 


2  k  ^  2cu-2k  ^ 

(1  -  xf)'^dx  <  C  Z  6  ^  / 

^  =0  -'d, 


2 

9  V , 


6/2<p<6  .  “1 

9x. 


xk  . 

x^)  dx 


k  2(r-2k  r  a  ^  ^ 

<C  Z  6  ^  L  L  a  ft  ft-  P 

1^1=0^  2 


a  2  -23^-2 ^.k.2 


sin  6  dpdedtp 


k  2a,-2k  r  2  -2y,+26,  -2My^2^-2 

<  C  Z  6  ^  I  /  I -I  r  P  ^  sin  ^  e 

0^=0  l£l=a^  6  a?3ra^ 


[4.7) 


2ii-i+k-2a,+2  2vu-2B^+2 

•  p  ^  sin  ^  6  sin  6  dpd6d(p 


<  C  Z  6 
01=0 


k  z'  ^  -2yi+2Bi  -2vi2'''2B2~2 

^  /  ft-  X  P  ^  ^  sin  ^ 

'i'=“4<p<6  3?3e23? 

-2a^  2vi,+2  2ii  +2 

•  (p  sin  6)  p  sin  6  sin  e  dpaBacp 


k  2a^-k-2oLj^+2  fnin(p^ ,  W2)+2 


<  C  6 


2  min(pj^,lJi2)+2 


Moreover , 


^  Iv  -  ^  ~f  ■  X5(P)Xg(P  sin  e))l^  dx 

P  V 

<cj  Ivl^  dx 

c  1  n 


“'p  sin  0<6 


(4.8) 


f.  -2u,  -2^2-2  2y,+2  2vu+2 

kcI  IvI  p  sin  e  P  sin  6  sin  Bdpdedcp 

•'p  sin  0<6 


z-  _  ~2yi  -2vi2+2  i  iMii  I V  Hi ,  M2 

<  c/  Ivr  p  ■‘■sin  6  (p  sin  6)  sin  G  dpdGdcp 

•/n  ci  rt 


2  niin(y^,y2)+2 


“P  sin  e<6 


<  C  6 


2  min(vij^,V^)+2 


For  t  >  0,  define 


K(v,  t)  .  inf  (llv  II  .  t  "''2"  k,,3-> 

vsVj^+V2  l-2(^  )  2  U  ; 


'  ■  r*  •-p 


For  0  <  t  <  1  let  ^  ~  ^  3*^*^  ^2  ~  '^6*  *^y  (^.7)  and  (4.8), 

min(y,  ,y_)+l  -k/2  +  min(y, ,  vi,)+l 

K(v,  t)  <C(6  ^  +  t  6  )  . 

?/k 

Choosing  6  =  t  ,  it  follows  that 


K(v,  t)  <C  t 


2  min(^^,l^)+2 
k 


0  <  t  <  1  . 


For  1  £  t  <  ®,  take  Vj^  =  v  and  V2  =  0  to  obtain 


K(v,  t)  <C  itvli  ,  , 

L2(I^) 


1  <  t  <  ® 


Hence , 


f  (r®K(v,  t))^^<cf 
*^0  ^  -70 


4min  (y,  ,^u)+4 
1  -26  +  - p:L_£ - 1 

t  dt 


+  C  llvll 


'  3/' 


t"  “  ^dt  <  » 


2min  (y,,y2)+2 

provided  that  0  <  6  < - -  .  By  the  definition  of  real  interpolation 

via  the  K-riiethod  together  with  [11,  Ineorem  2.1],  it  follows  that  for 
arbitrarily  small  e  >0, 

k  »  ^ 

2min  (y,,vu)+2-e  , 

=  z  ^  ^  (r) 


which  completes  the  proof  of  (i). 

The  proof  of  (ii)  is  essentially  the  same  as  that  of  part  (i)  except  that 


one  instead  sets  v^  =  vxx(p)- 


Having  established  the  regularity  of  the  solution  u  of  (A.l)  near  the 
edges  and  corners  of  ii,  the  piecewise  polynomial  approximability  of  u  in 
can  now  be  determined.  Let  1  £  £L,  oenote  the  corners 

of  it  and  let  0  <  i,  —  ^2  il  —  ***  eigenvalues  of  the 

Laplace-Beltrami  operator  on  with  homogeneous  Dirichlet  boundary 
conditions  on  where  is  the  portion  of  the  unit  sphere 
subtended  by  the  infinite  cone  which  coincides  with  fi  in  a  neighborhood  of 
o^.  Let  i  j  1.  denote  the  edges  of  and  let  u). 

denote  the  measure  of  the  interior  dihedral  angle  of  Ji  at  Yj* 
following  is  the  main  result  of  this  section. 

Theorem  4.1  Let  f  belong  to  C*(ii).  If  u  is  the  solution  of  (4.1)  then 
for  each  non-negative  integer  p  there  exists  e  P^(I2;A,D) 
such  that 

iiu  -  'll  It  1  <.C(u,  ^  ^ 

for  arbitrarily  small  e  >  0  where 

(4.9)  y  =  min 
1<JI< 
l<i<) 

and  C(u,  e)  is  independent  of  p. 

Pf ;  By  employing  an  appropriate  partition  of  unity,  the  approximability  of 
in  ft  may  be  determined  by  separately  considering  the  approximability  of  u 
near  the  edges  of  ft,  near  the  corners  of  ft,  ana  in  the  interior  of  ft. 

Since  u  is  infinitely  differentiable  in  the  interior  of  ft.  Theorem  2.2  may 
be  applied  with  any  value  of  s  yielding  an  arbitrarily  high  degree  of 
approximation  there.  Near  an  edge  y..  Lemma  4.2  and  Theorem  2.2  (with  an 
arbitrarily  large  value  of  s)  may  be  applied  to  the  terms  ^  and  h  of 


(4.4)  respectively;  this  together  with  Theorem  2.1  yields  an  approximation 

2  IT 

degree  of  -  e  near  Yj*  Near  a  corner  and  an  edge  Yj  with  as  an 
endpoint,  one  applies  Lemma  4.4  to  each  term  of  (4.5)  with  k'  =  [^],  k 
arbitrarily  large,  and  iJi2»  ^  chosen  as  follows: 


1.  For  the  term  p 


*  2  \  ^i,il  ■^4-' 


^m,n  '  take  =  -  2 


y2  =  1”  +  2n,  and  v  arbitrary  in  (i); 

Z  U)j 

2.  For  the  term  f^  take  y^^  =  k,  y2  =  ^  +  2n,  and  v  =  k'  in  (i); 

’1  ’  n - r  ‘'j 

~  2  '*’'V  fi.  4  1  I  1 

3.  For  the  term  p  '  ’  g^,  take  u  =  -  ■^  +  yj  a  *  t 

and  V  =  k  +  2  in  (ii); 

4.  For  the  term  h,  take  y  =  k  and  v  =  k'  in  (ii). 


Theorem  2.1  then  yields  an  approximation  degree  of  2 ‘min  | 

for  u  near  and  Yj*  The  proof  is  completed  by  adding  together  the 
various  approximating  piecewise  polynomials  obtained  in  this  manner  and 
applying  the  triangle  inequality. 

In  the  two-dimensional  result  Theorem  3.1,  the  degree  of  approximation  is 
dependent  solely  upon  the  measures  w.  of  the  interior  angles  at  the 
corners  of  the  domain  and  upon  the  boundary  conditions;  since  these  are  always 
known,  the  approximability  of  the  solution  of  (3.1)  is  readily  computed  in  all 
cases.  While  the  effect  of  singularities  along  the  edges  of  a  polygonal 
oomain  is  similarly  expressed  in  terms  of  the  dihedral  angle  measures  cu- 

«J 

in  Theorem  4.1,  the  presence  of  the  Laplace-Beltrami  eigenvalues  in  (4.9) 
makes  the  determination  of  the  approximation  degree  more  oifficult  in  the 
three-dimensional  case  since  the  relationship  between  these  eigenvalues  and 
known  quantities  in  the  problem  is  not  clear.  The  determination  of  the 
smallest  (or  fundamental)  eigenvalue  of  the  Laplace-Beltrami  operator 
on  various  domains  has  been  considered  by  a  number  of  authors  (see  e.g.  [8], 
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[20], [27]).  For  some  configurations,  can  be  obtained  analytically, 
however,  in  general  one  must  resort  to  numerical  techniques.  Figure  4.1  below 
depicts  a  solid  containing  a  notch  and  rectangular  crack.  Assuming  that 
(4.1)  is  posed  on  ii,  Table  4.1  lists  the  values  of  corresponding  to 
each  of  the  corners  of  0  together  with  the  method  used  to  compute  and 
the  quantity  -  +-^ ^  *  Considering  Table  4.1  together 

with  (4.9),  it  is  interesting  to  note  that  the  corners  of  type  A,  B,  or  C,  the 
degree  of  approximation  is  not  affected  by  the  "corner  singularities"  since 
these  are  weaker  than  the  singularities  along  the  edges  converging  at  such 
corners.  On  the  other  hand,  for  corners  of  type  D  or  E,  the  corner 
singularities  are  dominant. 


[ii 


[•j 


12.000 

Analytical  [27] 

% 

1.667 

Analytical  [8] 

>j  r 

S? 

1.500 

Analytical  [8] 

§  D 

0.453 

Numerical  [8] 

1 

0.596 

Numerical  [27] 

1 

Table  4.1 

5.  Numerical  Results 

The  relevance  of  any  asymptotic  error  analysis  such  as  the  one  presented 
here  depends  of  course  on  the  ability  of  the  underlying  method  to  achieve  the 
predicteo  asymptotic  order  within  practical  computational  limits.  Since  these 
issues  are  difficult  to  resolve  analytically,  one  must  rely  on  numerical 
experiments  to  determine  whether  it  is  possible  to  see  the  predicted 
asymptotic  behavior  of  the  error  without  using  an  unreasonable  number  of 
degrees  of  freedom.  Such  experiments  also  provioe  useful  information  about 
the  level  of  accuracy  that  can  be  expected  from  a  given  amount  of 
computational  effort.  In  this  section,  the  p-version  of  the  finite  element 
method  is  applied  to  two  model  problems  from  two-dimensional  linear  elasticity 
using  the  research  program  COMET-X  developed  at  the  Center  for  Computational 
Mechanics  at  Washington  University  in  St.  Louis  [7].  while  the  computations 
do  agree  with  the  asymptotic  orders  predicted  by  the  above  approximation 
theory,  the  main  conclusion  to  be  orawn  from  these  results  is  that  for 
problems  of  practical  interest,  the  p-version  does  in  fact  enter  the 
asymptotic  range  for  low  values  of  p.  Moreover,  it  is  seen  that  the  order  of 
convergence  is  completely  determined  by  point  singularities  in  the  solutions, 
and  that,  unlike  the  h-version,  the  point  at  which  the  p-version  enters  the 
asymptotic  range  is  virtually  unaffected  by  the  Poisson  ratio. 

In  order  to  link  the  following  computations  with  the  approximation  theory 
of  the  previous  sections,  a  number  of  preliminary  remarks  must  be  made. 

(i)  The  approximation  results  of  [11]  were  developeo  only  for  scalar 

problems.  However,  there  are  no  difficulties  in  generalizing  these 
results  to  the  approximation  of  vector-valued  functions  (such  as  the 
two-dimensional  displacement  fields  treated  below)  by  vectors  of 
piecewise  polynomials  in  the 


v< 


Sobolev  product  spaces  IH®  =  n  H®. 

Along  with  a  table  of  displacements  ana  stresses,  COMET-X  calculates 
the  strain  energy 


J(Up)  =  1/2  /jj  [C^^]'B[pyp]dfi 


of  the  computed  solution  ^  where  D  is  a  differential  operator 
and  B  is  an  elasticity  matrix  which  depends  on  the  elasticity 
modulus  £,  the  Poisson  ratio  v,  and  the  type  of  two-dimensional 
problem  being  considered,  i.e.  either  plane  strain  or  plane  stress 
(for  a  more  detailed  description  of  the  matrices  D  and  B,  the 
reader  is  referred  to  [28]).  In  terms  of  J,  one  defines  the  energy 


JJ-iJ  s  J(«) 


which  is  equivalent  to  the  usual  DH^  =  X  Sobolev  norm  on 
the  space  of  admissible  displacements  with  rigid  body  motions 
factored  out.  Since  the  error  in  the  energy  norm  corresponds  to  a 
least  squares  error  in  the  stresses,  it  is  generally  desired  in 
engineering  applications  to  achieve  5  to  10  percent  relative 
error  in  this  norm.  This,  of  course,  corresponds  to  a  .25  to  1 
percent  relative  error  in  energy. 

(iii)  In  two  dimensions,  (p+l)(p+2)/2  degrees  of  freedom  are  required  to 

uniquely  specify  a  polynomial  of  degree  p.  Since  the  p-version 

assumes  a  fixed  mesh,  it  follows  that  the  Dimension  N  of  the  finite 

2 

element  subspaces  grows  asymptotically  like  p  as  p  ®.  This 
implies  that  if  the  theory  of  the  previous  sections  preoicts  an  error 


.V  .V  .-V'  Ni 


estimate  of  the  form 


"  ii  ■  —  Cp“®,  a  >  0, 

where  u  is  the  exact  solution,  then  the  corresponding  estimate  of 
the  error  in  energy  with  respect  to  the  number  of  degrees  of  freedom 
is  given  by 

J(u  -  <  CN(p)"“. 

Due  to  a  well-known  orthogonality  property  of  the  error  u  - 
this  is  the  same  as 

(5.1)  I  J(u)  -  J(^)  I  <  CN(p)"“  . 

Consider  an  edge-cracked  square  panel  as  depicted  in  Figure  5.1.  In 
order  to  have  an  exact  energy  with  which  to  compare  energies  computed  by 
COMET-X,  let  it  be  assumed  that  the  tractions  i  =  1,  .  .  .,5  are 
such  that  if  u  is  the  restriction  to  the  panel  of  the  exact  solution  (as 
constructed  e.g.  in  [21, sec.  120])  for  an  infinite  cracked  panel,  then  u 
solves  the  finite  problem  of  Figure  5.1.  Near  the  crack  tip,  u  is  of  the 
form  r  »(0)  where  (r,e)  are  polar  coordinates  at  the  tip  of  the 
crack  and  £  is  a  smooth  vector  function.  Using  the  formula  for  ^  ,  one 
computes  J(u)  =  A2.16. 
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Figure  5.1:  An  edge-cracked  panel. 
Elasticity  modulus  E  =  1.0,  Poisson 
ratio  V  =  0.3,  Traction  o  as 
specified  in  the  text.  Plane  strain. 


Mesh  #  1 


Mesh  #  2 


Figure  5.2:  Two  meshes  for  the 
edge-cracked  panel  of  Figure  5.1. 
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Applying  the  p-version  to  this  problem,  one  subdivides  the  panel  in  some 
convenient  manner.  Two  possible  meshes  are  shown  in  Figure  5.2.  It  should  be 
noted  here  that  although  the  approximation  theory  presented  in  [11]  has 
assumed  triangular  subdivisions  this  is  not  essential,  and  the  results 
presented  are  equally  valid  for  quadrilateral  subdivisions  such  as  in  Figure 
5.2.  By  the  results  of  section  3,  one  obtains  that,  neglecting  e,  for  each 
polynomial  degree  p,  the  p-version  finite  element  solution  satisfies 


u  -  Mp  "ini  1  cp- 


where  C  is  independent  of  p.  Hence  the  error  in  energy  is  expected  to 
decrease  linearly  with  respect  to  the  number  of  degrees  of  freedom. 

Figure  5.3  shows  the  COMET-X  results  for  this  problem.  One  observes  that 
the  predicted  linear  convergence  is  achieved  and  that  the  asymptotic  range  is 
entered  for  p  as  low  as  4  for  either  mesh  #1  or  mesh  #2.  This  entry  into 
the  asymptotic  range  therefore  occurs  when  the  relative  energy  norm  error  is 
about  27  percent  and  14  percent  for  mesh  #1  and  mesh  #2,  respectively. 

The  superiority  of  mesh  #2  over  mesh  #1  is  clearly  due  to  the  fact  that,  given 
a  fixed  number  of  degrees  of  freedom,  mesh  y/2  allows  more  degrees  of  freedom 
to  participate  in  the  resolution  of  the  singularity  at  the  crack  tip  than  does 
mesh  #1.  Note  that  for  mesh  #2  a  relative  error  in  energy  of  less  than  one 
percent,  i.e.  a  relative  energy  norm  error  of  less  than  10  percent,  is 
achieved  with  p  =  6,  which  corresponds  to  482  degrees  of  freedom. 
Furthermore,  if  the  error  curve  for  mesh  #2  is  linearly  extended,  one  finds 
that  a  5  percent  relative  energy  norm  error  is  achievable  by  using  less 
than  2,000  degrees  of  freedom.  A  problem  of  this  size  could  be  handled  by  a 
large  mainframe. 
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Figure  5.3:  COMET-X  results  for 
the  edge-cracked  panel  of  Figure  5.1. 


Consider  next  the  edge-cracked  square  panel  in  Figure  5.4  where  the  lower 
side  of  the  crack  is  held  fixed  and  part  of  the  upper  edge  of  the  panel  is 
subjected  to  a  uniform  traction  as  shown.  Although  there  does  not  exist  an 
analytical  solution  to  this  problem,  it  can  be  rigorously  demonstrated  that 
near  the  crack  tip,  the  displacement  field  u  is  of  the  form 
where  (r,e)  are  polar  coordinates  based  at  the  crack  tip,  ^  is  a 
smooth  vector  function,  and  a  »  1/4.  Milder  singularities  of  a  similar 
form  exist  at  the  other  corners  of  the  panel.  By  the  approximation  theory  of 
section  3,  it  holds  that,  neglecting  c,  for  each  polynomial  degree  p,  the 
p-version  finite  element  solution  ul  satisfies 


Relative  error  in  energy 


V 


II  u  -  Up  Ij^l  < 


where  C  is  independent  of  p.  A  convergence  order  of  1/2  with  respect  to 
the  number  of  degrees  of  freedom  is  therefore  expected. 

Since  the  exact  energy  is  unknown  for  this  problem,  some  other  means  were 
required  to  obtain  a  reasonable  estimate  of  this  value.  For  this  purpose,  the 
self-adaptive  h-version  code  FEARS  (Finite  Element  Adaptive  Research  Solver) 
developed  at  the  University  of  Maryland  [5]  was  used.  Since  FEARS  provides  an 
aposteriori  estimate  of  the  error  which  may  be  employed  to  obtain  a  prediction 
of  the  exact  energy,  this  procedure  was  used  to  obtain  the  values  J(u)  =  1.57 
for  V  =  .3  and  J(u)  =  1.23  for  v  =  .49. 


t  t't  t 


Figure  5.4;  An  edge-cracked  panel 
with  constrained  crack.  Elasticity  modulus 
E  =  1.0,  Poisson  ratio  v  =  0.3  or  0.49, 
Traction  a  »  1.,  Plane  strain. 
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Figure  5.5  contains  the  COMET-X  and  FEARS  results  for  this  problem  using 
Poisson  ratios  of  v  =  .3  and  v  =  .49.  The  COMET-X  triangulation  of  the 
panel  is  shown  in  the  inset.  Since  FEARS  adaptively  selected  a  number  of 
meshes  to  obtain  the  results  of  Figure  5.5,  these  are  not  shown. 

One  observes  that  the  COMET-X  results  are  quite  similar  for  v  =  .3 
and  V  =  .49.  As  for  FEARS,  the  results  for  v  =  .49  are  somewhat  worse 
than  those  for  v  =  .3,  although  the  fact  that  an  asymptotic  rate  is  even 
observed  is  rather  significant  since  non-adaptive  h-version  codes  are  known  to 
perform  poorly  for  values  of  \>  near  the  incompressibility  limit  1/2.  The 
invariance  of  the  p-version  with  respect  to  the  Poisson  ratio  v  has  been 
previously  observed  by  other  investigators  and  is  rigorously  analyzed  in  [26]. 

The  convergence  order  for  COMET-X  is  1/2,  as  predicted,  whereas  the 
convergence  order  for  FEARS  is  1,  which  is  expected  due  to  the  fact  that  the 
adaptive  procedure  eliminates  the  effect  of  the  singularity  on  the  order  of 
convergence  leaving  only  a  dependence  on  the  polynomial  degree  (=1).  For  v 
=  .3,  the  COMET-X  error  is  less  than  that  of  FEARS  up  to  about  350  degrees 
of  freedom,  and  for  v  =  .49  the  COmeT-X  error  is  less  than  the  FEARS  error 
up  to  about  2500  degrees  of  freedom.  These  cross-over  points  correspond  to 
a  relative  energy  error  of  about  9.5  percent  and  3.7  percent, 
respectively.  Thus,  in  very  low  accuracy  ranges,  the  COMET-X  results  are 
somewhat  better  especially  for  the  Poisson  ratio  v  =  .49.  On  the  other 
hand,  in  the  higher  accuracy  ranges  FEARS  achiieves  a  lower  error  than 
COMET-X.  In  fact,  by  linearly  extending  the  error  curves  of  Figure  5.4,  one 
easily  checks  that  in  order  to  obtain  results  which  have  a  relative  error  of 
10  percent  in  the  energy  norm,  if  v  =  .3  then  COMET-X  would  require 
about  35,000  degrees  of  freedom  as  opposed  to  3,000  for  FEARS.  If  v  = 

.49  then  COMET-X  would  require  about  25,000  degrees  of  freedom  to  about 


Comat-X: 


4  5  6  7  8 


v  *  0.49 


V  a  0.49 


Number  of  degrees  of  freedom 


Figure  5.5;  COIET-X  and  FEARS  results 
for  the  edge-cracked  panel  of  Figure  5. A. 


7,000  for  FEARS.  Both  of  the  latter  numbers  of  degrees  of  freedom  for 
COMET-X  must  be  regarded  as  being  outside  practical  computational  limits. 


In  summary,  although  the  above  numerical  results  are  in  no  way 
exhaustive,  they  do  suggest  the  following: 

(i)  The  order  of  convergence  of  the  p-version  agrees  with  that  predicted  by 
the  approximation  theory.  Moreover,  this  order  of  convergence  is  completely 
controlled  by  a  principal  singularity  occurring  at  a  point (s)  on  the  boundary 


(ii)  The  asynptotic  range  is  achieved  for  low  values  of  p.  Furthermore,  the 
asymptotic  range  has  generally  been  entered  by  the  time  the  relative  error  in 
the  energy  norm  has  reached  10  to  30  percent. 

(iii)  The  p-version  is  not  sensitive  to  changes  in  the  Poisson  ratio  v, 
even  when  v  is  close  to  the  incompressibility  limit. 

(iv)  If  the  dominant  singularity  is  of  the  form  r°^  with  a  _>  1/2  , 

then  it  is  possible  to  achieve  a  relative  error  in  the  energy  norm  of  5 
percent  without  exceeding  practical  computational  limits.  If,  however,  a 
<  1/2,  say  e.g.  a  =  1/4  as  in  the  second  example,  then  it  may  not  be 
possible  to  obtain  a  relative  energy  norm  error  of  less  than  20  percent 
without  using  an  excessive  number  of  degrees  of  freedom.  For  problems  with 
such  strong  singularities,  it  seems  clear  that  an  adaptive  p-version  or  a 
combination  of  h-  and  p-versions  [3]  will  be  necessary  to  obtain  the  level 
of  accuracy  required  for  many  problems. 

Of  course,  the  above  considerations  have  avoided  other  aspects  which 
figure  significantly  in  the  actual  cost  of  finite  element  analysis.  These 
include  the  sparsity  of  the  resulting  matrices,  the  complexity  of  the  input 
data,  issues  of  data  structure  and  data  management,  etc.  Questions  such  as 
these  are  discussed  in  [6]  and  [22-24]. 
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